#ifndef UHM_MAT_LINAL_OPS_HXX_
#define UHM_MAT_LINAL_OPS_HXX_

#define UHM_MAT_LINAL_OPS_UNROLL 4

namespace uhm {

  // Mat Ops
  // -------
  void random_matrix(Mat_Linal_ &A);
  void random_spd_matrix(Mat_Linal_ &A);

  void identity_matrix(Mat_Linal_ &A);
  void norm1(Mat_Linal_ &A, Double_ &val);
  void norm_inf(Mat_Linal_ &A, Double_ &val);
  
  void copy(Mat_Linal_ &A, Mat_Linal_ &B);
  void copyt(Mat_Linal_ &A, Mat_Linal_ &B);
  void axpy(Int_ alpha, Mat_Linal_ &A, Mat_Linal_ &B);



  // LINAL Interface
  // ---------------
  Int_ apply_diag_scal     ( Int_ side, Int_ conj, Int_ uplo,
                             Mat_Linal_ &A );

  Int_ apply_diag_scal     ( Int_ side, Int_ conj,
                             Mat_Linal_ &A, Mat_Linal_ &B );
  Int_ apply_diag_inv_scal ( Int_ side, Int_ conj,
                             Mat_Linal_ &A, Mat_Linal_ &B );

  Int_ apply_pivots( Int_ side, Int_ trans, Mat_Linal_ &p, Mat_Linal_ &A );

  Int_ chol_nopiv  ( Int_ uplo, Mat_Linal_ &A );
  Int_ chol_piv    ( Int_ uplo, Mat_Linal_ &A, Mat_Linal_ &p );
  Int_ chol_incpiv ( Int_ uplo, Mat_Linal_ &A, Mat_Linal_ &p );

  Int_ ldl_nopiv   ( Int_ uplo, Mat_Linal_ &A );
  Int_ ldl_piv     ( Int_ uplo, Mat_Linal_ &A, Mat_Linal_ &p );
  Int_ ldl_incpiv  ( Int_ uplo, Mat_Linal_ &A, Mat_Linal_ &p );

  Int_ lu_nopiv    ( Mat_Linal_ &A );
  Int_ lu_incpiv   ( Mat_Linal_ &A, Mat_Linal_ &p );
  Int_ lu_piv      ( Mat_Linal_ &A, Mat_Linal_ &p );

  Int_ trmm  ( Int_ side, Int_ uplo, Int_ trans, Int_ diag,
               Int_ alpha, Mat_Linal_ &A, Mat_Linal_ &B );
  
  Int_ trsm  ( Int_ side, Int_ uplo, Int_ trans,
               Int_ diag, Int_ alpha, Mat_Linal_ &A, Mat_Linal_ &B );
  Int_ gemm  ( Int_ transa, Int_ transb, Int_ alpha,
               Mat_Linal_ &A, Mat_Linal_ &B, Int_ beta, Mat_Linal_ &C );
  Int_ gemm_l( Int_ transa, Int_ transb, Int_ alpha,
               Mat_Linal_ &A, Mat_Linal_ &B, Int_ beta, Mat_Linal_ &C );
  Int_ gemm_u( Int_ transa, Int_ transb, Int_ alpha,
               Mat_Linal_ &A, Mat_Linal_ &B, Int_ beta, Mat_Linal_ &C );
  Int_ syrk  ( Int_ uplo, Int_ trans, Int_ alpha,
               Mat_Linal_ &A, Int_ beta,  Mat_Linal_ &C );
  Int_ herk  ( Int_ uplo, Int_ trans, Int_ alpha,
               Mat_Linal_ &A, Int_ beta,  Mat_Linal_ &C );

  // UHM Ops
  // -------
  void random_matrix(UHM_Linal_ &U);
  void random_spd_matrix(UHM_Linal_ &U);

  void identity_matrix(UHM_Linal_ &U);
  void schur_x2r(UHM_Linal_ &U);

  void copy(UHM_Linal_ &UA, UHM_Linal_ &UB);
  void copy_A(UHM_Linal_ &UA, UHM_Linal_ &UB);
  void copy_B(UHM_Linal_ &UA, UHM_Linal_ &UB);

  void axpy(Int_ alpha, UHM_Linal_ &UA, UHM_Linal_ &UB);
  void axpy_A(Int_ alpha, UHM_Linal_ &UA, UHM_Linal_ &UB);
  void axpy_B(Int_ alpha, UHM_Linal_ &UA, UHM_Linal_ &UB);
  
}


#endif
